NUMERICAL EVALUATION OF THE OSCILLATORY 
INTEGRAL OVER cxp{iTrx)x^/'' BETWEEN 1 AND INFINITY 



RICHARD J. MATHAR 



Abstract. Real and imaginary part of the limit 2N — >■ oo of the integral 
J-^ exp{iTTx) y/xdx are evaluated to 20 digits with brute force methods after 
multiple partial integration, or combining a standard Simpson integration over 
the first half wave with series acceleration techniques for the alternating series 
co-phased to each of its points. The integrand is of the logarithmic kind; 
its branch cut limits the performance of integration techniques that rely on 
smooth higher order derivatives. 



1. Scope 

1.1. M. R. Burns' Constant. 

Definition 1. The MRB constant is the sum of the series 'lE', A037077] 

2N 

m 

N- 



(1) M= lim y(-l)"V^ = y(-l)'=(fci/'=-l)« 0.18785964. 

n=\ k=l 

Direct summation of the alternating series is slow and generates roughly 3 valid 
digits after ten thousand terms (Table [1]). Euler summation [1] (3.6.27)] [10] is 
successful in accelerating the convergence, witnessed in Tabled] 

Table 1 . Partial sums of ([1]) as a function of the upper limit of summation. 

10 0.313231759254. 

10^ 0.211329543346. 

10^ 0.191323989712. 

10-* 0.188320351076. 

10^ 0.187917210140. 



More efficient methods lead to even quicker convergence, as demonstrated in 
Table [3] An accuracy of 60 digits is reached after 100 terms and will be sufficient 
for all purposes of this script. 
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Table 2. Approximations to ([T]) after Euler resummation of the 
first k terms. 

k Ejl) 

10 0.f878858861f380073035f 382438464680824292327407645248188116946. 
20 0.187859649854050194658445181421685949965109596411421113134589. 
40 0.187859642462068655529781630996634559217485607603915205816966. 

100 0.187859642462067120248517934054273314215151463271236583869269. 

200 0.187859642462067120248517934054273230055903094900138786171982. 



Table 3. Approximations to ([T]) with the first Cohen- Villegas- 
Zagier algorithm using terms up to fc [3]. 

k M 

10 0.187859642389333316567457476113016727048599369932998191180459. . . 

20 0.187859642462067119674255755542940758484982176117045969528027. . . 

40 0.187859642462067120248517934054273140023454509840554949525330. . . 

100 0.187859642462067120248517934054273230055903094900138786171986. . . 



1.2. Oscillatory Integral. The integrated analog of the series is a complex- valued 
integral of oscillatory character, which is difficult to evaluate by direct integration 
if the upper limit becomes large, illustrated by Figure [T] 

Definition 2. (Sequence of oscillatory integrals) 

Not convergent in the continuum limit at TV — >■ cxd, the limit of the sequence of 
integrals with an integral difference in the upper limits 2N exists. The objective of 
this work is to evaluate this limit Af/. 

Definition 3. (Ultraviolet limit of the sequence) 
(3) Mi = lim I{2N). 

N~>-oo 

Remark 1. The absolute value |Af/| w 0.6876523689 is close to M + ^ [H 
A157852]. Changing the upper limit to 2N + 1 increases Mj by 2i/7r. 

To compute M/, the manuscript looks at repeated partial integration to quench 
the integrand at large x in preparation for standard methods of sampling along 
the abscissa (Sections 12.11 and I2.2p . investigates splitting the integral into an alter- 
nating series and a base interval (Section I2.3p . expansion of x^/^ into a series over 
(log x/x)" fSection l^^ . changing the path of integration in the complex plane (Sec- 
tion [^3]), and considers reverse application of the Euler-Maclaurin integral formula 
(Appendix [C|. 

2. Numerical Analysis 



2.1. Iterated Partial Integration. A partial integration of ([2]) yields 

j-2N 

(4) / e^'^'^a^^/'^da; = --e^^x^/^ 



I / 1 1 - log X 



27V 

TT , TT 




The limit iV — >■ oo can be performed in the pre-integrated term, 
(5) M, = + 1 r e-^'/'^l^^dx, 

which essentially compresses the oscillations with a factor cx log(a;)/a;^, as shown 
in Figure [5J 
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Figure 2. Real and imaginary part of the integrand in ([5]). 
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A second partial integration of JS]) may follow, 
(6) f"" e'--x^/-l^:}^dx 



^...^i/. l-log^ ^ f + 2{x - 1) log(x) + V o: ^^ 



with the limit 



Repeating, a simple scheme for the n-th derivative of the base function /, 
can be phrased as a set of coefficients a, 

log^(x) 

r,s>0 



(9) /(")(a;)=:ri/- ^ a„,,^. 



Explicit computation with the chain rule establishes the recurrence 

(10) an+l.r,s = an,r-2,s - Ol.n,r-2,s-l + (s + l)an,r-l,s+l - {t " l)a„,r-l,s- 

The initial conditions are 

(11) ctn,r.,s =0 if r < or s < 0; ao,o,o = 1- 

Equations ([5]), ([7]) and the representation through rt-fold partial integration are 
summarized with 

(12) Mi^b{n)+i^\ j e''''' f^''\x)dx, 
defining pre-integrated terms 

(13) h{n + l) = h(n)--f^^-'\l); 6(0) = 0; 6(1) = --. 

The infinite interval from 1 to oo can be mapped onto the interval from to 1/2 
with the substitution (one of a family of rational maps) 



1 1 , 

1 + a; y 



1 /2 

(15) / e""/(")(x)dx= / e*""(^)/^"Ha;(y))^ 



y2 



A combination of and (ITSI) yields: 



Algorithm 1. Compute Mj with a Simpson integration on s points in the interval 

0<y< 1/2.- 

(16) Mj = bin)+(^^^ ^'^%--(J/)/(«)(a;(y))^. 
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Table 4. Algorithm [TJ Simpson integration with s points in the 
interval < j/ < 1/2 and the n-th partial integration (fT2|) inserted 
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32000 


07077603931155557088 


-0 68400038943794081422 


5 


64000 


0.07077603931152496616 


-0.68400038943793426656 


6 


2000 


0.07077603931115051881 


-0.68400038943952383435 


6 


4000 


0.07077603931156321151 


-0.68400038943794556862 


6 


8000 


0.07077603931154254337 


-0.68400038943792386069 


6 


16000 


0.07077603931152730519 


-0.68400038943793206498 


6 


32000 


0.07077603931152869782 


-0.68400038943793231267 


6 


64000 


0.07077603931152878997 


-0.68400038943793211168 


8 


2000 


0.07077603932853563807 


-0.68400038944555242777 


8 


4000 


0.07077603931259171964 


-0.68400038943840850746 


9 


2000 


0.07077603925904942234 


-0.68400038951437536286 


9 


4000 


0.07077603930824878310 


-0.68400038944270983422 



Table |4] illustrates that a choice of n near 5 or 6 yields optimum convergence (be- 
cause h{n) then approximate Mi best), and that integration with s « 60 000 abscissa 
points generates of the order of 13 valid digits. Note that Romberg (Richardson) 
extrapolation with the standard 15 : 1 weighting of step width halving does not 
work as the integrand is not in the polynomial class. 

Remark 2. In a variant of this mapping on a finite interval, the transformation 
X — 1/u in the n-th partial integration yields 

(17) Mj - b{n) + {i/nr C e^-/"/(")(l/u)^. 

Jo " 

The precision is worse than with Algorithm]^ by approximately one digit. I have not 
looked into advanced schemes for this type of oscillatory integrals [71 [TT] or Sidi 's 
generalized methods of extrapolation. 

2.2. Exponential Scaling. A characteristic of Algorithm[T]is that the integrand is 
basically reduced by another factor 1 /x for each additional partial integration. The 
variable transformation logx — z, x — e^ , dx = e^dz, helps to achieve exponential 
scaling as the integration variable heads towards infinity, at the cost of an irregular 
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chirp factor in the complex exponential: 

(18) Fm= f e''^^a;i/^ ^~^Qg'^ rf^^ / (l _ z)e^^°'^P(^)e^<='^P(-^)-^dz. 



Based on 

(19) y gZ7roxp(z)+2^^ ^ _lgi7rexp(z)^ 

and "borrowing" a factor in ^TE\\ in the integrand, 

POO 

(20) (l-z)e*''°'^P(")+^e^<=''P(-^)~2z^^^ 

^ log m 

a partial integration generates 
(21) 

= '-tn .^i^i^ + i r g.-cxp(.)g.cxp(-.)-2. _ 3 + _ 2z + z2 

Alternatively, this resuhs applying the substitution logx = z to (IT2|) . After this 
scheme of partial integrations has been repeated n times, the non-oscihating expo- 
nential factor in the integrand is exp[2; exp(— z) — (n + l)z]. 

Algorithm 2. Perform n partial integrations of il8\) . then use another transfor- 
mation u ~ e^^, z = — logw to map the range 0<z<ootoO<u<lin the 
remaining integral, and integrate this over u with a Simpson method. Eventually 
insert this Fi in 

(22) Mi = -- + -F^ 

TT TT 



as s 



een in l^)- 



An accuracy of 10^^^ can be reached sampling two million points (Tabled and 
is reported in the summary. Comparison with Table 0] shows that one to two digits 
are gained relative to Algorithm [1] 

2.3. Longman's Method. An integral F with an undulating trigonometric factor 
multiplied by a monotonous g{x) may be split into an integral over the half period 
with an alternating series attached to each point in that interval [l4l [T3l [5] . 

(23) F,^ = / e"^5(x)dx e*'^(™+'+'').g(m + I + y)dy 

= (-)" e^^y ^(-)'ff(m + I + y)dy. 

The Z-series is only alternating if the function g{x) is monotonous and does not 
change sign. In addition, Euler resummation assumes that the series converges. 
With M/, x^/^ has a maximum at a; = e, so we integrate over 1 < a; < 3 with any 
other method, setting m = 3, and relay by one partial integration, g{x) = f^^^x), 
to feature a g{x) that has a single sign with decreasing |s('(a;)| in [m, oo). 
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Table 5 . Convergence of Algorithm [2] integration of the n-th 
partial integration of Fi with s equidistant points in the interval 
< u < 1. 



n 


s 


3fi(M/) 


9(M/) 


2 


4000 


0.07077612979610666804 


-0.68400038256040301228 


2 


8000 


0.07077604264870771610 


-0.68400037471801246748 


3 


4000 


0.07077603954212465077 


-0.68400039095784542350 
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-0 68400038943798882953 
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K 
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5 


16000 


07077603931151256105 


-0 68400038943794236597 
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6 


4000 


0.07077603931156163745 


-0.68400038943790335638 


6 


8000 
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-0.6840003894379321291827445 
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-0.68400038943794252308 


9 


4000 
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Algorithm 3. Compute Mj by Filon-Simpson-integration of Fm in i23\) over the 
interval < y < 1 fAppendixlA]) and 
(24) 



Mr 



-F 



g{x) = xi/^ 



1 — log X 



(to 



With this choice, g{x) has a maximum near x ~ 4.3, associated with the 
zero of f'^^^x), which (after detailed inspection) does not destroy the alternat- 
ing property — if higher order /^"^(a;) were employed, to would have to be chosen 
differently. 

The speed of convergence of the method is demonstrated with Table IH] As the 
number of evaluations of g is the product of n and I, it turns out effectively to be 
slower than Algorithm [2j 



2.4. Taylor Series the Logarithmic Term. The most advanced method expands 
the non-oscillatory term of ^ into the Taylor series of the exponential: 
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Table 6. Convergence of Algorithm |3] with a Simpson integration 
on n abscissa points over [0,1], truncating the alternating series 
after the l-th term followed by extrapolation [3]. 
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32 


70 


0.07077603721021819390 


-0.68400038753980049654 


64 


70 


0.07077603918043104863 


-0.68400038931936816305 
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70 


0.07077603930333884842 


-0.68400038943052296458 


256 


70 


0.07077603931101698843 


-0.68400038943746907333 


512 


70 


0.07077603931149681599 


-0.68400038943790318846 


1024 


70 


0.07077603931152680433 


-0.68400038943793032039 


2048 


70 


0.07077603931152867859 


-0.68400038943793201613 


4096 


70 


0.07077603931152879573 


-0.68400038943793212212 


8192 


50 


0.07077603931152880305 


-0.68400038943793212874 


8192 


60 


0.07077603931152880305 


-0.68400038943793212874 


16384 


50 


0.07077603931152880351 


-0.68400038943793212916 


16384 


60 


0.07077603931152880351 


-0.68400038943793212916 



Table 7. Convergence of the partial sum of Algorithm H) 

maxn jR(Mj) Q(Afj) 

1 0.05762490298863188764 -0.68331060191932132015 

2 0.06935454902524610824 -0.68451362283943263006 

3 0.07066781734932533318 -0.68408393964817446557 

4 0.07076978736326279015 -0.68400839361204835470 

5 0.07077575770475264223 -0.68400096184084805332 

6 0.07077602957520227166 -0.68400042266805241181 

7 0.07077603908225272182 -0.68400039107408096452 

8 0.07077603931050590528 -0.68400038950815469850 

9 0.07077603931177817100 -0.68400038944060977421 

Algorithm 4. (Logarithmic expansion o] x^^^) 

n2N p2N p2N 1 1 n 

(25) / e^'^^xi/^^dx = / e"^e^>°s^dx= / e"^[l + V -7-^^^]da; 
Ji Ji Ji -, n'. 



--+1:-^/ 



n>l 
2N 



log" X 



The integrals are the V{TT,n,n) defined in equation pip in Appendix [BI Sum- 
ming over values taken from Tables IMTT] produces Table [71 An accuracy of 10~^^ 
in real and imaginary part of Mj requires summation up to n = 15 — an estimation 
derived from results of Remark |4] — , and has not been worked out. 

2.5. Contour Deformation. The path of the integration may be deformed to a 
straight line (hypotenuse) from z = 1 to z = 2iV(l -I- ri) with adjustable slope 
r > towards the real axis, and a straight line (short leg) parallel to the imaginary 
axis back to 2A^ on the real axis. The contribution of 5z to expiinz) leads to a 
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Table 8. Algorithm [5] Simpson integration along the line from 
z = 1 to z = 2A^(1 + ri) with s abscissa points. 
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exponential reduction of the integrand as the distance to the real axis grows. (More 
complicated paths appear not to be more efficient.) The short leg of this triangle 
contributes 

(26) lim / e*'^^+'°s^/2(i2 = 

to the integral. 
Algorithm 5. Compute 

■ .2N{l+ri) 

(27) Mi = + lim / e''''+^°s'/'dz 

with a Simpson integration on s points Zj — I + j(l + ri) At, j = 0, . . . s. 

The dependence on three configuration parameters makes quality assessment 
more dilficult than with the other methods, illustrated by Table[Sl The contribution 
missing for any finite A'^ is estimated by 

(28) / e"''+^°«'/'dz^ e'^^'dz^ -e^'''^^^+'^\ 

J2N{1+Ti) J2N{1+Ti) 

which serves as a guideline how large r ought be made given a targeted accuracy 
and an upper limit N . 
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3. Summary 

The value of Mi is 

(29) lim / e^'^^xi/^dx 

w 0.0707760393115288035395- 0.68400038943793212918i. 

The integral features highly oscillatory behavior and a logarithmic factor in the 
main integrand, which sets up an interesting test case outside the range of methods 
that assume "nice" analytic properties in the complex plane. 



Appendix A. Filon-Simpson Rule 

The Simpson rule of integration of a function G{y) is an interpolation between 
three abscissa points (j/o, G(0)), (yi, G(l)) and (1/2, G(2)) by a quadratic polynomial 
iny [3 (25.2.1)], 

^ (y-yi)(y-2/2) ^.„^ , (2/-2;o)(2/-y2) , (y-yo)(y-yi) 

[vo - yi)\yo - y2) (yi - yojivi - y2j {y2 - yajm - yi) 

followed by integration of this polynomial in the limits ya < y < y2- A refinement 
in our case, based on the same quadratic interpolation, is 

e'''yGm{y)dy = TTrr^h'^ + Z/i^tt^ - Sihir + (2 - /i7ri)e2"'*]G(0) 

- ^^[-* + h7T + ii + h7T)e^'^'']Gil) 

- :Tr^[2 + + (-2 + 2h^T:^ + 3i/i7r)e2'"'']G(2), 

supposing equidistant abscissa points yi — yo + h and 2/2=2/0 + 2/i. This explicit 
recognition of the exponential factor gains roughly one additional digit in accuracy 
in Table [6] compared to the evaluation with e"^ incorporated in the value of G. 



Appendix B. Fichtenholz Integrals 

B.l. Fundamental Form. In this appendix, a set of integrals V{a, k, s) is targeted 
as an aid to Algorithm |4l This is basically evaluating the generalized integro- 
exponential function [161 112] for a complex-valued parameter z = —ia. 

Definition 4. (Generalized Integra- Exponential Function) 
(31) V{a,k,s)= / e'^'^-^dx. 



Remark 3. There are three simple extensions: 

• Cases with a scaling factor b in the logarithm can be reduced to this funda- 
mental form by binomial expansion of {log b + log x)'^ , 

(32) f e-^^dx = g log'=-'(6) F(a, I, s). 
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Reading 



(33) / e-'^^dx = b^-^ re-Z^'^dy, 



y 

(obtained through the substitution bx ^ y) from right to left shows that 
other lower limits than 1 are also accessible once the V are known for 
general a. 

• Integer powers of sines or cosines at the place of the exponential lead back 
to the fundamental form via Euler 's formula: 

sm"\ax)^dx = ^ ^ (-l)'F[a(™ - 21), k, s]; 



(35) 



cos"^{ax)^dx = — Y,r]v[aim-2l),k, 

^ 1=0 ^ 



Real and imaginary part of the value V{a, 1, 1) 
(36) re'--}2i^dx^ r coa{ax)^^dx + I H sm{ax)^-^dx 



are computed separately. The imaginary part is 

(37) rsHax)'-^dx^ rMa.)'^dx- [\Hax)'-^dx 



with one constituent ^ (4.421. 1)][1 (865.63)][S] 
(38) 



f°° lo2 X n 

/ sin(ax)-^da; = (7 + log a) w -2.704825746060380848849568 (a = tt). 

Jo X 2 

The integral over [0, 1] is evaluated by Taylor expansion of the sine which leads to 
the well converging representation 

(39) / sin(ax)^d2; = -a^(-l)«— — — -— — ^ 
Jo X ^ (2?i + l)!(l + 2n)^ 

= -a2F3 3/2^3/2^3/2 I ^ -2-6581349165086 (a = n). 

The difference between this value and (|38)) represents ((37|) . 



(40, si„(ax)^dx . --h + log„) + „i:(-')- (,„^i),,i + ,„). 

n>0 

The value at a = vr is the head entry in Table ID 
The real part is started from [9, (3.761.9)] 

(41) rSq(^rf.= ^Meos(^), 
which is differentiated with respect to ji with the product rule, 

(42) / cos(ax) ^ dx = ^^^^ i'lbiu) cos — — In a cos — — — sin —\ . 
. n x^-f af V ' 2 2 2 2 7 
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Table 9. Table of the imaginary part of V{tt, 1, s). 

_s 5y(7r, l,s) 

1 -0.046690829551739977074516092264 

2 -0.050400599397438879223041567776 

3 -0.044723677797644192936988882199 

4 -0.036181141258573216997609321919 

5 -0.027931519676734642467423612590 

6 -0.021096986691682143229642825070 



Table 10. Table of the real part of V{Tr, 1, s). 

_s jRF(7r, l,s) 

1 0.057624902988631887643485542240 

2 0.029913203983934978439301792236 

3 0.010937363639874260291206201403 

4 -0.000250069139610211209861368961 

5 -0.006024230915536561502482189260 

6 -0.008508918812024751462009533761 



The reflection formulas for the r-function and Digamma function ijj show that in 
the limit fi^O [I, (6. 1.17), (6.3. 7)] 

(43) cos(ax)-A_dx ^ ^ + Hi+l"") + ^ + ^^^^ 

where 7 « 0.57721 is the Euler-Mascheroni constant [TFl A001620]. Expansion of 
the cosine in the standard Taylor series [HI (1.411.3)] and interchange of integration 
and summation yields 

(44) r^rf-=E 



Differentiation with respect to [i, introduces the logarithm, 
(45) / cos(aa:) i_^ dx = — y 



n>0 



(2n)!(2n + /i)2 



We subtract this from (^5)) and notice that the ~ 1/ //^-singularity cancels with the 
term n = of the series as /i — >■ 0, 



(46) cosiax)^dx = -—+-fi^j+\naj+^ + Y^ 



logx , 7r2 /T , \ In^a (-a^)" 



(27i)!(2n)2 

The sum converges quickly, the value at a = tt is the first entry in Table 1101 
Combined with pp]) this reads 

log a; , TT^ /7 1 \ In^ a tt . , , , v-^ (iaY' 



(47) / e-^i^f:d. = -^+7Q+lna)+i^-|.(7 + loga) + 5: 



l«2 
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which constitutes a "root" value V{a, 1, 1) in the tree of integrals ([3T|) . A summary 

of (|42l) and the associated imaginary part is 

(48) 

rt!(n + 



e-'^d. ^ Via, 1, 1-,) ^ ^Me^W^ (^(,) - In a + f ) +^ 

^ ^ Tl>0 



B.2. Higher Powers of the Rational. Integration of (|T7)) with respect to a is a 
measure to increase the parameter s, the power of x in the denominator, by one: 

(49) r da' r d.e--'^ = H - l]'-^d. = H e^^'-^d. - L 
Jo Ji X Ji Wi X * 

(50) 



,.,,.loga; , . 



9 1 2 

-— + 7(^ + Ina - 1) + ^ - Ina + 1 - y (7 + Ina - 1) 



^ (n + l)!n2' 

n>l 



where [T7| 



(^a)"+^- _ (a»)^+^- / 1,1,1 I 



n>l 

Inserting a — tt yields the second lines in Table [10] and [9l The concept of (|49| 

generalizes to higher powers k of the logarithm, 

(52) 

da'Via', fc, s) = - / [6"^"= - 1]— I— = -iy(a, fc, 1 + s) + iT/(0, fc, 1 + s). 
« Ji 2:^+'* 

The last term, the non-oscillatory integrals at a = 0, are known 9, (4.272.6)], 

(53) F(0,0,s) = ^, mfc,^)- 

Milgram's equation (2.29) [16]. Iterated application of this rule computes the chain 
of V{a, 1, 1) ^ V{a, 1, 2) V{a, 1, 3) ^ ... as follows: 

(54) ^ ^ ^ 

r e'^'^^^dx = 4? - ^(70 + 2a In a - 3a) - ^(In^ a - 31na + 7/2) + i 
Ji a;3 48 4 4 4 

3, . v-^ (m)"+2 



+ ^«(7 + ln-:^)+a*+2^ 



(55) 



(n + 2)!n2' 

n>l ^ ^ 



e-^d^ ^ __^_l_^(^+2 lna-ll/3)+y ^(--+- In a-- In^ a)+-z 



TTa\ , , , a2 1 ^ (ia)"+3 

-(7 + lna-ll/6)- — + - + ^ 



12 " ' ' 2 9^ (n + 3)!n2' 

The cases of a = tt are in Tables [9] and [10] 
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B.3. Higher Powers of the Logarithm. As seen in Section TB.H differentiation 
with respect to the parameter /i increases the power of the logarithm: 

1 fc+i 
,v,^m 



(56) 



ds 



V{a, k,s) — — 



-dx 



-V{a,k + l,s). 



To support differentiation with respect to the s-parameter, we recompute V{a, 1, 2- 
fj), which we re-integrate over a as in ([42^ : 



(57) da cos(a x) dx = / sin(Qa:) ^_ dx 

Jo Jo ^ ^ Jo X ^^ 

r(M) 



a''-i(l - ^) 



1-M 



in a cos — — sm — — 



The integration is applied in parallel to the complementary interval < a; < 1, 
sin(ax)— — -dx — 



(58) 







X 



E 

n>0 



(2n+ l)!(2n + M)2' 



and the difference is 
(59) 

r(^) 



svii{ax)—^^dx = '^V{a, 1, 2 — /i) 



a^-i(l - ^) 



x^ 



1 



1-M 



— In a COS sm — 



n>0 



(2n+ l)!(2n + M) 



2 ■ 



This is differentiated with respect to /i, and the singularities ^ 2a//i'^ from the 
incomplete Gamma-function and the n-sum cancel in the limit /it — >■ 0: 



(60) 



sin(ax) 



log^ X 



dx — a 



7(-7 + 2 - In a) In a + Y^7r2(ln a - 1 + 7) - -((3) 



1 T 

21na + In^ a In^ a - ^ 27 + 7^ 

3 3 



The explicit value at a — tt is 
(61) 

loe^ X 

sm{Trx)-^dx = 5F(7r, 2, 2) 



-2«E 



n>l 



(2n+ l)!(2n)3' 



-0.00240604184022261982751961704408. 



This demonstrates the technique. Starting from V{a, 2, 2), a ladder of integrals is 
constructed according to ([52l) . Each time, a term ^ {—l)'^kl{iay~^/^''~^^ cancels — 
carried over from the simple pole of the F-function ()41[) through k differentiations 
and s — 1 integrations — when the complementary integrals of < a; < 00 and 
< X < 1 are combined. Numerical examples are gathered in Table [TT] 

Remark 4. Partial integration 



(62) 



. log'' X , f log*" X 1 



-dx 



dx 



Tog' 



fc+i 



, log'' X 



1 fc-i 
log X 



log'' X 



log xdx 
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5y(7r,fc,s) 



2 2 0.0234592920732284411947929 

2 3 0.0147167653246107850628908 

2 4 0.0080788243950624846223234 

2 5 0.0038282314868382609838401 

2 6 0.0013846381369360916659413 

2 7 0.0000998710536555974796592 

2 8 -0.0005097889780015996357674 



-0.0024060418402226198275196 
-0.0078739079798083225909317 
-0.0087094002295813942921939 
-0.0076206659960747444416239 
-0.0060334458936144149013900 
-0.0045501919692512628703906 
-0.0033548439064072548162349 



3 3 0.0078796099444753496396155 

3 4 0.0053790452810071493354663 

3 5 0.0032308011441634551812315 

3 6 0.0017674846795167938450865 

3 7 0.0008825074832995147577407 

3 8 0.0003872374479338897353530 



0.0025780991475489869676813 
-0.0004578919913424486140616 
-0.0014908072125213319989936 
-0.0015900953082259670812147 
-0.0013496431705418864940445 
-0.0010417218132785702486827 



4 4 0.0024472803344989672143380 

4 5 0.0018067919115537410157224 

4 6 0.0011430051029650516946321 

4 7 0.0006599851739729340692718 

4 8 0.0003564954008069045751359 



0.0018131048670266608466124 
0.0004312822675128498086469 
-0.0001368354074002352164914 
-0.0003012860412794823958833 
-0.0002990689541546756726797 



5 5 0.0007164409787822497618537 

5 6 0.0005826627556204031506092 

5 7 0.0003859348621737006202841 

5 8 0.0002304000284892623408819 



0.0008918125440361658853376 
0.0003133177357343641402474 
0.0000540098043633385948340 
-0.0000404135204291844446708 



6 6 0.0001957467237331888336746 
6 7 0.0001826715069173610144948 
6 8 0.0001272475588547481432200 



0.0003882044128618852155461 
0.0001565025492374147392079 
0.0000473426202252359704144 



yields a contiguous relation for three points in a triangle in the square grid of {k, s)- 
pairs (Milgram's equation (2.4) jlfij ); 

in Q — 1 

(63) V{a, k, s) = —^V{a, k + l,s-l) + ^^^(a, ^ + 1, «)■ 

Since the calculation of values at large k is a laborious task, this formula offers a 
route to cheaper calculation by (i) tabulation ofV{a,k,s) at some small k up to a 
rather large s, involving only integrations of powers oflna multiplied by powers of 
a [9l (2.722)], (ii) numerical calculation of the V{a,k + Ak, s) up to the desired k 
with some brute force method like the \/x mapping {11^ , which converges well since 
s is large, (Hi) telescoping from s backwards with ^63) to fill the table for increasing 
k and decreasing s. 

Remark 5. The last term in \6S\) can be eliminated via 156]) . 

(64) (1 + k)V{a, k, s) = -iaV{a, + 1, s - 1) + (1 - s)-^V"(a, fc, s), 

ds 

The solution of this inhomogeneous differential equation is 

(65) V{a, k, s) = \y+k (^~*" J k + l,s~ l)(s - l)''ds + const^ . 
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Appendix C. Inverse Euler-Maclaurin 

A standard idea of integration is to split the integral over intervals that are 
commensurable with the frequency of the oscillation, to replace the generic factor 
g in the integral by some approximation which allows integration in closed form — 
assuming that a massive cancellation can be obtained — , then to gather the sum 
over the intervals with some Euler-Maclaurin approach [2 . Applied to ([23| . 

I^ca oo fk+l 

(66) F™= / e'^'-'g{x)dx= ^ / e'^-g{x)dx, 

g is approximated by its Taylor series, [151 18]. 

(67) g(^)=^_gW(fc)(^_fc)<i. 

oo „i oo -J 

(68) F„= Yl e"'=y_ d2/e^-^^-5('^)(fc)/. 

Definition 5, (Moments of Filon Quadratures) 

(69) S{d)^ j'^e^^yy'^dy- S{d) ^ d = 0, 1, 2, 3, . . . 

Initial values are 

(70) 5(0) = 0; 3(1) = -. 

TT 

The recurrence is 

i 



(71) Sid + I) = -{[! + i-l)'] + id+l)Sid)} 



TT 

With 9j (3. 761. 5), (3. 761. 10)] the cases for even and odd indices are 



(72) Si2d+l) = 2zl^^--^(-.^)^ 



3=0 



(2j + 1) 



(73) Si2d) = 2UlY^^_L_(^^,y_ 



OO oo CXD 



This rewrites 
(74) 

F™= E e*"'=E5(d).g('^)(fc) = (-l)'"+i5]5]5(%('')(m+2/+l). 

k=m+l,m+3,m+5,... d=0 1=0 d=l 

The Euler-Maclaurin formula proposes to replace the sum over the d-th derivatives 

by 

(75) E.g('*)(m + 2/ + l) = - .g('*)(m + l)+ / g^'^\x)dx 

y _?^2^+V'+''Hm + l)V 



15=1,3,5. 
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Fm becomes a double sum over d and D, and rcsummation proposes Algorithm [6] 
to accumulate the derivatives of the (smooth) function g to calculate 

i , ^ 1/^1 — log X 

(76) Mi^-iF-2); g = x'/^ . 

Algorithm 6. (1-sided Fourier- Euler-Maclaurin) 

2(-)'"+iF = -5(l)5(o)(m + l) 

oo / V(d-l)/2\ \ 

+ Y^[s{d)-S{d+l)~ Y: S{d-l-2l)^^2^+^^9^^Hm+l). 

Implementation of this approach reveals that the sum over the d-th derivatives 
shows converging behavior only up to d « 6. I attribute this to the same loga- 
rithmic branch cut that constraints the useful depths of the partial integrations in 
Algorithms [D and El 
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